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Abstract 

We propose a semi-discrete finite difference multiscale scheme for a concrete corrosion model con- 
sisting of a system of two-scale reaction-diffusion equations coupled with an ode. We prove energy and 
regularity estimates and use them to get the necessary compactness of the approximation estimates. Fi- 
nally, we illustrate numerically the behavior of the two-scale finite difference approximation of the weak 
solution. 
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1 Introduction 

Biogenic sulfide corrosion of concrete is a bacterially mediated process of forming hydrogen sulfide gas and 
the subsequent conversion to sulfuric acid that attacks concrete and steel within wastewater environments. 
The hydrogen sulfide gas is oxidized in the presence of moisture to form sulfuric acid that attacks the 
matrix of concrete. The effect of sulfuric acid on concrete and steel surfaces exposed to severe wastewater 
environments (like sewer pipes) is devastating, and is always associated with high maintenance costs. 

The process can be briefly described as follows: Fresh domestic sewage entering a wastewater collection 
system contains large amounts of sulfates that, in the absence of dissolved oxygen and nitrates, are reduced 
by bacteria. Such bacteria identified primarily from the anaerobic species Desulfovibrio lead to the fast 
formation of hydrogen sulfide (H2S) via a complex pathway of biochemical reactions. Once the gaseous 
H2S diffuses into the headspace environment above the wastewater, a sulfur oxidizing bacteria - primarily 
Thiobacillus aerobic bacteria - metabolizes the H2S gas and oxidize it to sulphuric acid. It is worth noting 
that Thiobacillus colonizes on pipe crowns above the waterline inside the sewage system. This oxidizing 
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process prefers to take place where there is sufficiently high local temperature, enough productions of hy- 
drogen sulfide gas, high relative humidity, and atmospheric oxygen; see section 2.2 for more details on the 
involved chemistry and transport mechanisms. Good overviews of the civil engineering literature on the 
chemical aggression with acids of cement-based materials [focussing on sulfate ingress] can be found in 
[3, 12, 17, 22]. 

If we decouple the mechanical corrosion part (leading to cracking and respective spalling of the con- 
crete matrix) from the reaction-diffusion-flow part, and look only to the later one, the mathematical problem 
reduces to solving a partly dissipative reaction-diffusion system posed in heterogeneous domains. Now, 
assuming further that the concrete sample is perfectly covered by a locally-periodic repeated regular mi- 
crostructure, averaged and two-scale reaction-diffusion systems modeling this corrosion processes can be 
derived; that is precisely what we have done in [9] (formal asymptotics for the locally-periodic case) and 
[21] (rigorous asymptotics via two-scale convergence for the periodic case). 

Here, our attention focusses on the two-scale corrosion model. Besides performing the averaging proce- 
dure and ensuring the well-posedness of the resulting model(s), we are interested in simulating numerically 
the influence of the microstructural effects on observable (macroscopic) quantities. We refer the reader to 
[5], where we performed numerical simulations of such a two-scale model. Now, is the right moment to 
raise the main question of this paper: 

Is the two-scale finite difference scheme used in [5] convergent, i.e., does it approximate the weak solution 

to the two-scale system? 

It is worth mentioning that there is a wealth of multiscale numerical techniques that could (in principle) 
be used to tackle RD systems of the type treated here. We mention at this point three approaches only: 

(i) the multiscale FEM method developed by Babuska and predecessors (see the book [8] for more Refs.), 

(ii) computations on two-scale FEM spaces [16] / two-scale Galerkin approximations [19, 18], and (iii) the 
philosophy of heterogeneous multiscale methods (HMM) [7]. We choose to employ here multiscale finite 
differences (multiscale FD) mimicking the [two-scale] tensorial structure present in (ii). Our hope is to 
become able to marry at a later stage the two-scale Galerkin approximation ideas from [19, 15] in a HMM 
framework eventually based on finite differences. Our standard reference for FD-HMM idea is [1]. 

The paper is structured as follows: Section 2 introduces the reader to the physico-chemical background 
of the corrosion process, two-scale geometry, and setting of the equations. The numerical scheme together 
with basic ingredients like discrete operators, discrete Green formulae, discrete trace inequalities etc. are 
presented in section 3, while the approximation estimates together with the interpolation (extension) and 
compactness steps are the subject of Section 5 . We conclude the paper with Section 6 containing numerical 
illustrations of the discrete approximation of the weak solution. 

2 Background and statement of the problem 

2.1 Two-scale geometry 

We consider the evolution of a chemical corrosion process (sulfate attack) taking place in one-dimensional 
macroscopic region Q. := (0,L), L > 0, that represents a concrete sample along a line perpendicular to the 
pipe surface with x = being a point at the inner surface in contact with sewer atmosphere and x = L being 
a point inside the concrete wall. Since we do not take into account bulging of the inner surface due to the 
growth of soft gypsum structures, the shape of the domain Q. does not change w.r.t. the time variable t. 

We denote the typical microstructure (or standard cell [11]) by Y := (0,1), £> 0. Usually cells in concrete 
contain a stationary water film, and air and solid fractions in different ratios depending on the local porosity. 
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Generally, we expect that, due to the randomness of the pores distribution in concrete, the choice of the 
microstructure essentially depends on the macroscopic position x G il, i.e., we would then have Y x ; see [20] 
for averaging issues of double porosity media involving locally periodic ways of distributing microstructures, 
and [9] for more comments directly related to the sulfatation problem where pores are distributed in a locally 
periodic fashion. In this paper, we restrict to the case when the medium Q. is made of the same microstructure 
Y periodically repeated to pave perfectly the region. Furthermore, since at the microscopic level the involved 
reaction and diffusion processes take place in the pore water, we choose to denote by Y only the wet part of 
the pore. Efficient direct computations (with controlled accuracy and known convergence rates) of scenarios 
involving Y x as well as the corresponding error analysis are generally open problems in the field of multiscale 
numerical simulation. 

2.2 Chemistry 

Sewage is rich in sulphur-containing materials and normally it is without any action on concrete. Under 
suitable conditions like increased temperature or lower flow velocity oxygen in sewage can become depleted. 
Aerobic, purifying bacteria become inactive while anaerobic bacteria that live in slime layers at the bottom 
of the sewer pipe proliferate. They obtain needed oxygen by reducing sulfur compounds. Sulfur reacts with 
hydrogen and forms hydrogen sulfide (H2S), which then diffuses in sewage and enters sewer atmosphere. 
It moves in the air space of the pipe and goes up towards the pipe wall. Gaseous H2S (further denoted 
as H2S(g)) enters into the concrete pores (microstructures) via both air and water parts. H2S(g) diffuses 
quickly through the air-filled part of the porous structure over macroscopic distances, while it dissolves in 
the thin, stationary water film of much smaller, microscopic thickness that clings on the surface of the fabric. 

There are many chemical reactions taking place in the porous microstructure of sewer pipes which de- 
grade the performance of the pipe structure depending on the intensity of the interaction between the chemi- 
cal reactions and the local environment. Here we focus our attention on the following few relevant chemical 
reactions: 



Dissolved hydrogen sulfide (further denoted as H2S(aq)) undergoes oxidation by aerobic bacteria living in 
these films and sulfuric acid H2SO4 is produced (reaction (la)). This aggressive acid reacts with calcium car- 
bonate (i.e., with our concrete sample) and a soft gypsum layer (CaS04 ■ 2H2O) consisting of solid particles 
(unreacted cement, aggregate), pore air and moisture is formed (reaction (Id)). 

The model considered in this paper pays special attention to the following aspects: 

(i) exchange of H2S from water to the air phase and vice versa (reaction (lc); 

(ii) production of gypsum at micro solid-water interfaces (reaction (Id)). 

The transfer of H2S is modeled by means of (deviations from) the Henry's law, while the production of 
gypsum is incorporated in a non-standard non-linear reaction rate, here denoted as 77; see (6) for a precise 
choice. Equation (7) indicates the linearity of the Henry's law structure we have in mind. The standard 
reference for modeling gas-liquid reactions at stationary interfaces (including a derivation via first principles 
of the Henry's law) is [6]. 



H 2 S(aq) + 20 2 -> 2H+ + SO 2 " 
10H+ + SO4 2 + organic matter ->• H 2 S (aq) + 4H 2 + oxidized matter 

H 2 S(aq)^H2S(g) 
2H 2 Q + H+ + SO 2 .- + CaCQ 3 -> HCO3- + CaSQ 4 • 2H 2 Q. 



(la) 
(lb) 
(lc) 
(Id) 
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2.3 Setting of the equations 



Let S := (0,r) (with T £ (0,°°)) be the time interval during which we consider the process and let £2 and 
Y as described in section 2.1. We look for the unknown functions (mass concentrations of active chemical 
species) 

u\ : £2 x S — > R - concentration of H2S(g), 

«2 : £2 x T x S — > R - concentration of H2S(aq), 

1(3:^x7x5^1 - concentration of H2SO4, 

U4 : £2 x S — >• R - concentration of gypsum, 

that satisfy the following two-scale system composed of three weakly coupled PDEs and one ODE 

d t u\ — d\d xx u\ = d2d y U2\y=o, in £2, (2a) 

d t U2 —d%dyyii2 — — £(u2,u-j), in £2 x Y, (2b) 

d t UT, —d-idyyUi = £ (112, 113), in £2 x Y, (2c) 

<9,m 4 = 77 (m 3 m 4 ), in£2, (2d) 

together with boundary conditions 

u\ = uf, on {x = 0} x S, (3a) 

d\d x ui — 0, on {x = L} x S, (3b) 

-d 2 dyU2=B^(Hui-U2), on£lx{y = 0}xS, (3c) 

t/ 2 <3y«2 =0, on £2 x {y = £} x S, (3d) 

— d-idyii^ ~Q, on £2 x {y = 0} x S, (3e) 

didyiii = — 7](m3, M4), on £2 x {y = x 5, (3f) 



and initial conditions 



Mi =«5, on£2x{f = 0}, 

M2 = t<2, on £2 x F x {r = 0}, 

m 3 = M3, on £2 x Y x {t = 0}, 

1/4 = W4, on£2x{f = 0}. 



(4) 



Here t4, k £ {1,2,3}, are the diffusion coefficients, Bi is a dimensionless Biot number, H is the Henry's 
constant, a,j3 are air-water mass transfer functions, while rj(-) is a surface chemical reaction. Note that 
Ui (i = 1, ... ,4) are mass concentrations. Furthermore, all unknown functions, data and parameters carry 
dimensions. 



2.3.1 Technical assumptions 

The initial and boundary data, the parameters as well as the involved chemical reaction rate are assumed to 
satisfy the following requirements: 

(Al) d k > 0, k £ {1,2,3}, Bi M > 0, H > 0, uf > are constants; 
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(A2) The function £ represents the biological oxidation volume reaction between the hydrogen sulfide 
and sulfuric acid and is defined by 

C:K 2 ^M, CM :=ar- /3s, (5) 

where a,j3 e 

(A3) We assume the reaction rate r] : R 2 — > R + takes the form 

. fkR(r)Q(s), forallr>0,s>0, 

7]{r,s) = < (6) 
10, otherwise, 

where k > is the corresponding reaction constant. We assume 77 to be (globally) Lipschitz in both argu- 
ments. Furthermore, R is taken to be sublinear (i.e., R(r) < r for all r £ R, in the spirit of [2]), while Q 
is bounded from above by a threshold c > 0. Furthermore, let R € W 1,oo (0,Mi) and Q S W l '°°(0,M4) be 
monotone functions (with R strictly increasing), where the constants M3 and M4 are the L°° bounds on M3 
and, respectively, on 114. Note that [5, Lemma 2] gives the constants M^,M4 explicitly. 

(A4) u\ e// 2 (n)nL7(n), e [L 2 (n;// 1 (y))] 2 x x r)] 2 , u\ G// 1 (n)nL7(n); 

2.3.2 Micro-macro transmission 

Terms like 

Bi M (Hu l {x,t)-~u 2 (x,y = 0,t)) (7) 

are usually referred to in the mathematical literature as production terms by Henry's or Raoult's law; see [4]. 
The special feature of our scenario is that the term (7) bridges two distinct spatial scales: one macro with 
x £ £1 and one micro with y G Y. We call this micro-macro transmission condition. 

It is important to note that in the subsequent analysis we can replace (7) by a more general nonlinear 
relationship 

In that case assumption (A2) needs to be replaced, for instance, by (A2') 

0$ e C 1 ([0,Mi] x [0,M 2 ];R), 38 globally Lipschitz in both arguments , (8) 

where M\ and Mi are sufficiently large positive constants 1 . Note that a derivation of the precise structure of 
S3 by taking into account (eventually by averaging of) the underlying microstructure information is still an 
open problem. 

2.4 Weak formulation 

As a next step, we first reformulate our problem (2), (3), (4) in an equivalent formulation that is more suitable 
for numerical treatment. We introduce the substitution ii\ := u\ — uf to obtain 

d t ii\ — d\d xx Ui = d2dyU2\ y= Q 1 in £2, (9a) 

d t u 2 - d2d yy ii2 = -C("2,«3), in£2xT, (9b) 

d t uj —d$d yy Ui = £(u2,U3), in £2 x Y, (9c) 

d t U4 = Tj(w3|y=^W4), in£2, (9d) 



Typical choices for Mi ,Mi are the L°°-estimates on u\ and H2; cf. [5] (Lemma 2) such M\,M2 do exist. 
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together with boundary conditions 



u x = 0, on {x = 0} x S, (10a) 

did x ti\=0, on {x = L}xS, (10b) 

-d 2 d y u 2 =Bi M (H(u l +u^)-u 2 ), on £2 x {y = 0} x S, (10c) 

d 2 d y u 2 = 0, onQ.x{y = £}xS, (lOd) 

-d 3 d y u 3 =0, on£2x{y = 0}xS, (lOe) 

d^dyUi — — T}(«3,M4), on £2 x {y = ^} x 5, (lOf) 



and initial conditions 



u\ = u\ — uf =: fij, on x {f = 0}, 

mt = j<2) on n x 7 x {f = 0}, 

M3 = z<3, on£2xTx{r = 0}, 

M4 = «4, on£2x{f = 0}. 



(11) 



We refer to the system (9), (10), (11) as problem (P). Also, for the ease of notation, we denote u\ again as 
u\ and vt\ as u\. 

Now, we can introduce our concept of weak solution. 

Definition 1 (Concept of weak solution). The vector of functions (u\, u 2 ,u 3 ,un) with 

mi EL 2 (S,H^Q.)), (12) 

d t meL 2 (SxQ.), (13) 

Ui eL 2 (S,L 2 (a,H l (Y))), ie{2,3}, (14) 

d tUi eL 2 (SxQ.xY), ie{2,3}, (15) 

u 4 (-,x,y) £H l (S) 7 for a.e. (x,y) ESlxY, (16) 

is called a weak solution to problem (P) if the identities 



and 



d,ui<pi+di / d x uid x (pi = / <3 y M 2 |v=o<Pi, 
> Jn Ja 

I I d,U 2 (p2+d 2 / / d y U 2 dy<p 2 = - / / C("2)"3)<P2- / <9 y «2|v=0<P2, 

iiiiy JnJy JnJr Ja 

/ d t U3(p3+d 3 d y U3d y (p3= / C(m 2 ,M 3 )^3 - / J7(H3lv=^"4)<P3, 

<9(M 4 = Tj(M3| y= f,M4), 



hold for a.e. t E S and for all <p := (<pi , 92, <Ps) G ff,} (fl) x [L 2 (n;i/' (7))] 
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We refer the reader to [5, Theorem 3] for statements regarding the global existence and uniqueness of 
such weak solutions to problem (P); see also [19] for the analysis on a closely related problem. 
The main question we are dealing with here is: 

How to approximate the weak solution in an easy and efficient way, consistent with the structure of the 
model and the regularity of the data and parameters indicated in (A1)-(A4)7 
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3 Numerical scheme 



In order to solve numerically our multiscale system (2)-(4), we use a semi-discrete approach leaving the time 
variable continuous and discretizing both space variables x and y by finite differences on rectangular grids. 
In the following paragraphs we introduce the necessary notation, the scheme and discrete scalar products 
and norms. 

3.1 Grids and grid functions 

For spatial discretization, we subdivide the domain £2 into N x equidistant subintervals, the domain Y into N y 
equidistant subintervals and we denote by h x := L/N x , h y := £/N y , the corresponding spatial step sizes. We 
denote by /; the vector (h x ,h y ) with length \h\. 



be, respectively, grid of all nodes in £2, grid of nodes in £2 without the node at x = (where Dirichlet bound- 
ary condition will be imposed), grid of all nodes in Y, grid of nodes located in the middle of subintervals of 
£lh, and grid of nodes located in the middle of subintervals of Y/,. Finally, we define grids ©/, := £2/, x Yf, and 



Next, we introduce grid functions defined on the grids just described. Let Sf/, := {m/, | Uf, : £1), —> R}, 

:= {«/, | m/, : £2^ — > R}and S% := {v/Jv/, : Q. e h — > R} be sets of grid functions approximating macro variables 
on £2. Let J^), := {«/, | «/, : CO/, — ► R} and Jff, := {v/ ? | v/, : (sf h — > R} be sets of grid functions approximating 
micro variables on £2 x Y. These grid functions can be identified with vectors in R , whose elements are the 
values of the grid function at the nodes of the respective grid. Hence, addition of functions and multiplication 
of a function by a scalar are defined as for vectors. 

For m/, € Sf/, we denote m '.= Uh(x\), and for «/, S J£), we will denote utj '■= Uk(xi,yj). For v/, S §h we will 
denote v, +I/2 := v/,(x /+ i/ 2 ), and for v A G we will denote v,- j+l / 2 := v/, (jc i -,y >+1/2 )- 

We frequently use functions from restricted to the sets £2/, x {y = 0} or £2/, x {y = £}. For m/, G 
we will denote these restrictions as m/,| v= o and w/, | v= £, and we will interpret them as functions from i.e., 
u/,\ y =o G % and M A | y= ^ e g},. 

3.2 Discrete operators 

In this section, we define difference operators defined on linear spaces of grid functions in such a way they 
mimic properties of the corresponding differential operators and, together with the scalar products defined 
in Sec. 3.4, fulfill similar integral identities. 

The discrete gradient operators V/, and V v /, are defined as 



Let 




{xi:=ih x | i = 0,...,N x }, 
{ Xi | i=l,...,N x }, 

!.v ( : A / W\. 

{x i+1/2 :=(i+l/2)h x \i = 0,,..,N x -l}, 
{y J+ i/2:=U+l/2)h y \i = 0,...,N y -l}, 



©| := £l h x Y e h . 



V h : % -). 4, 



(Vv/,M A ). ■ 




u h G & h 
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while the discrete divergence operators div/, and div v /, is 



div A :<f ft ^°, 
div v/i : J^-)- &h: 



(div A v A )j := 
(div y /,v A )y := 



v. . , 1 — v. . 1 



V/, G 



The discrete Laplacian operators A/, and A y i, are defined as A/, := div/, V/, : §^ — > <&? and A,,/, := div v /, V y /, 
J£/, — >• J£},, i.e., the following standard 3-point stencils are obtained: 



(A h u h ] 



(AyhUh)ij 



-2«,- + «,- +1 
A? 

"ij-l + Uj.j+ 1 



To complete the definition of the discrete divergence and Laplacian operators, we need to specify values 
of grid functions on auxiliary nodes that fall outside their corresponding grid. At a later point, we obtain 
these values from the discretization of boundary conditions by centered differences. 



3.3 Semi-discrete scheme 

We can now construct a semi-discrete scheme for problem (2). Note that we omit the explicit dependence 
on t and we interchangeably use the notation and ii/, for denoting the derivative of «/, with respect to t. 



Definition 2. A quadruple {u\ 1 uj v uf v uf l } with 



ul4&C\[QJ];%) and i? h ,*? h zC l ($,T\\& h ) 



2 „3 n \ , 



is called semi-discrete solution of (2), if it satisfies the following system of ordinary differential equations 



du\ 
dt 

dt 
dw| 

dt 
dut 

dt 



d x A h u\ - Bi M {H{u\ + u\ } ) - ul\ y=0 ) , 
d 3 Ay h ul + C(u 2 h ,ul), 



r\{ul\ y= i,u\), 



on flj, 
on ft)/,, 
on ft)/,, 
on i2/,, 



(17a) 
(17b) 
(17c) 
(17d) 
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together with the discrete boundary conditions (i = 0, . . . ,N X ) 

4=0, (18a) 

di\((yhu\) Nx ^+{V h ul) Nx _^ = 0, (18b) 

-di\ (p?jk4\-\ + ( V v/>4), i ) = BP(H{u}+u?) - ulo), (18c) 

rf 2 i ((VtfHj&^+i + (V,*^) w _i) = 0, (18d) 

((Vy,4). _ i + i) = 0, (18e) 



(18f) 



and the initial conditions 

4(o) = ^u?, 4(0) = ^, 
4(o) = ^ 2 4 4(o) = ^, 

where 5^ and are suitable projection operators from i2 to i2/ ; and from Q. x F to ft)/,, respectively. 



(19) 



Remark 1. The boundary conditions (18b)-(18f) are a centered-difference approximation of conditions (3) 
and are written so as to stress the relation between the two. Using the definition of discrete V/ ; and V T /, 

follows: 



operators, (18) can be rewritten in terms of auxiliary values of uL k= 1, ... ,4, on nodes outside the grids as 



4 


= o, 




U N x +\ 


= 4*- 


-1 • 


4-i 


=4i 


+ ^lBi M (H(uj+u? 

"2 


l i,Ny+\ 


= u Ih 


,-l> 


3 

4-i 


3 

— Mj-j. 




'Lvy + l 


= u In : 


,-i--^7(4iV"i)- 



(20a) 
(20b) 

&)> (20c) 

(20d) 
(20e) 

(20f) 

Proposition 3. Assume (A1)-(A4) to be fulfilled. Then there exists a unique semi-discrete solution 

{4,4,4,4} e C\[0,T];& h ) x c\[0,T];^ h ) xC\[0,T];& h ) xC\[0,T];& h ) 
in the sense of Definition 2. 

Proof. The proof, based on the standard ode argument, follows in a straightforward manner. □ 

3.4 Discrete scalar products and norms 

Next, we introduce scalar products and norms on the spaces of grid functions ^/,, §h, J^j,, Sf/, and we show 
some basic integral identities for the difference operators. 
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Let (yD?= and (7j)jL be such that 



1 l<i<N x -l, 



1 l<j<N y -l, 



* : \i /e{0,Ak}, ' ^ \± je{OA}, 

and define the following discrete L 2 scalar products and the corresponding discrete L 2 norms 



(uf,,Vhh h '-=h x Y Yi"iVi, 

Xi<E£l h 

\\ u h\\% : = J {u h ,Uh)y h , 



U h ,V h GW h , 
u h G %, 



(21) 

(22) 
(23) 



u h ,v h ew%, 
u h e Stf , 



(24) 
(25) 



(u h ,v h )j? h :=h x h y Y YiYjUijVij, 
xtjecoh 



Uh,V h 6 



(26) 
(27) 



(u h: v h )g h :=h x u /+i/2Vi+i/ 2 , 

Xi+l j 2 ea e h 



IKIk ■= \J{u h ,\x h )g h 



u/,,v/, e S"h, 



(28) 

(29) 



{^h,Vh)jtr h -=h x hy Y ^ u ij+i/2 v /J+i/2! 



H||^, := ^/(u/,,u A )^, 



U/, € J*£/j. 



(30) 



(31) 



It can be shown that a discrete equivalent of Green's formula holds for these scalar products as well as 
other identities as is stated in the following lemmas. 



Lemma 4 (Discrete macro Green-like formula). Let Uf, € @h and V/, £ S'h such that 

no = 0, un x +i — un x -\, 
v N x +i/2 = —Vn x -i/2- 



Then the following identity holds: 



(32) 
(33) 



(34) 
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Lemma 5 (Discrete micro-macro Green-like formula). Let Uh G &h and V/, G ^7, 5mc/i that 

-\ (vfc.-i^ + v^i/a) = 3t, i (^.-1/2 + ^+1/2) = 5 A 2 , (35) 
«* ,-i = + 2/i y 5/ , Uk,Ny+i = Uk,N y -i + 2h y 8%, (36) 
for i = 0,... ,N X , and Si , 8% G Sf/,. 77ie« the following identity holds: 

(u h ,diY yh v h )^ h = -(VyhUh,Vh)jg> h + {uh\y=Q,8l)% + {i*h\y=Ny,8l)y h . (37) 
We also frequently make use of the following discrete trace inequality: 

Lemma 6 (Discrete trace inequality). For «/, G ^ f/iere exists a positive constant C depending only on £2 
such that 

'/.!>•=* Ik < C(||u A ||.^ + IIV^MftH^). (38) 



Proof. Our proof follows the line of thought of [10]. We have that for Uf, G 

Ny-l Ny 

|"<yv>,| < 52 - «y I + £ ^l"y|- 

Squaring both sides of the inequality, we get 

(u,, Ny ) 2 <Ai + Bi, (39) 

where 

/N y -\ \ 2 /Ny \ 2 

A ' : = 2 ( L I M '' J+ 1 — "oil and B i '■ = 2 ( E ^ I M '7 1 J • 

Applying the Cauchy-Schwarz inequality to A,-, we obtain 

Ny-l / X 2Ny-l Ny-l / \ 2 

A,- < 2 £ ,„ (^±pA I A, = 2£ £ (^^) • 

7=0 \ / 7=0 ;=0 V "> / 

Similarly, using the Cauchy-Schwarz inequality we get for Bj 

Ny Ny Ny 

B, < 2 £ Ay(uy) 2 £ T 2 /^ = 2£ £ ^(u,y) 2 . 
.7=0 ;=0 7=0 

Multiplying (39) by ^J 1 /^, summing over / G {0, . . . ,N X } and then using the bounds on A, and Bj, it yields 
that: 

N x ( N x N y -l 2 \ 

Y,yl h M,Ny) 2 <2t\ £ £ ^w, + £ r!7}h x hy(uij) 2 ), 

i=0 \('=0 7=0 '2/ *y6<ty, / 

that is 

1 1 «/,!>■=£ \\% h < c (llVvAM/,11 2 ^ + \\u h \\%^j , 

from which the claim of the Lemma follows directly. □ 
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4 Approximation estimates 



The aim of this section is to derive a priori estimates on the semi-discrete solution. Based on weak 
convergence-type arguments, the estimates will ensure, at least up to subsequences, a (weakly) convergent 
way to reconstruct the weak solution to problem (P). 



4.1 A priori estimates 

This is the place where we use the tools developed in section 3. 

In subsequent paragraphs, we refer to the following relations: From scalar product of (17a) with q>l £ 
(17b) and (17c) with <p ; 2 e ^7, and <p%, respectively, and (17d) with <p£ £ % to obtain 

(ul<Phhg=M A hul<Phhg -Bi M {Hul-u 2 h \ y= o,(ph)&o, (40) 

("ft. <Ph)* h = dl{&yhul,<pl)& h - a(u 2 h , (?l), 9h +P(ul, (p 2 h ),? h , (41) 

("i 9ftV* = <fe (Ayhulvl)^ + a(ul (pl)^ -p(ul (pl)jr h , (42) 

("t9ftk = (n( u h\y=lA),9h) % - (43) 

Note that ul and Vhtyl satisfy the assumptions of Lemma 4, u\ and Vy,9? satisfy the assumptions of Lemma 
5 with 5/ = ^-{Hu\ - u\ Q ) and <5 A 2 = 0, and u\ and V v/! ^ with 5/ = and 5 A 2 = -^("j^y Thus, 
using Lemmas 4, 5 and properties of the discrete scalar products we get 

{u\,(pl) % +d l (V h u 1 h ,V h ( P 1 h )^ = -Bi M (Hul-ul\ y= o,<pl) rfii , (44) 

("ft, (ph)# h + d2(Vyhuly y h<Ph)M = Bi M (Hu\ - 4\ y=0 , (p 2 \ y= oh h -a(u 2 h , (p 2 h ),? h + j5 {u\,(pl)& h , (45) 

(ul Vh)j? h + d^ yh uly yh yl), m = - (77 {u\\ y=l , 4), 9ft \y=l)cg h + «("ft> 9ftK ~ P ( u l 9ftK > (46) 

(4>9ftk = {n{ul\y=eA),<Ph) % - (47) 

Lemma 7 (Discrete energy estimates). Let {u\ 1 iP' h ^ir' jv u' h \ be a semi-discrete solution of (2) for some T > 0. 
77ien /f holds that 

max (|K(0|||, + \\u 2 h (t)\& h + H(t)\\% h + H(f)\\k) < C > (48) 
r 



l|Vft^||^ + PywlW^ + \\V yh ul\\^ h )6t < C, (49) 

where C := C(j|w^(0) ||^ + ||m|(0)||^ + ||«^(0)||^ + II m a(0)II^, *w'f/i C £>e/ng a positive constant indepen- 
dent ofh x , ky. 

Proof In (44)-(47), taking (9 A \ 9^, 9^, 9s) = {u\,u\,u\,ufy, summing the equalities, applying Young's 
inequality on terms with (m?,m?)^-., dropping the negative terms on the right-hand side, and multiplying the 
resulting inequality by 2 give 



d 
d7 



(Kill, + ll"ftll^, + IK'llk + Kill,) +2rfi||V A 4||| > + 2rf 2 ||V^||2 ri+2rf3 ||v >ftM 3||2 
< -2Bi M {Hu\ -u 2 \ y=0 ,ul) % +2Bi M (Hul- U 2 \ y=0 ,u 2 \ y=Q ) % 

+ Ci\\u 2 h \\ \ + C, 1 1 1 1 \ + 2(x](u\ \ y=( , u 4 h ), u 4 h - u\ \ y=e ) „ , 
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where C\ := a + j3 > 0. Expanding the first two terms on the right-hand side we get 

-2{Hu\ - u 2 h \y = a,u\)c<? h +2(Hu\ - u\\ y= ^u\\ y= ^) % = -2H\\u\\\ 2 % 

+ 2(1 +H){ul,ul\ y=0 ) % -2\\u 2 h \ y= 4l h < + ((1 +H)e-2)\\ul\ y=0 \\l h , 

where we used Young's inequality with e > 0. Choosing e sufficiently small, the coefficient in front of the 
last term is negative, so we have that 

-2(Hu h -ul\ y= o,u h )& h +2(Hu h -ul\ y= o,ul\y=o)% < Cz\\u h \\\, 

where C 2 ■= > 0, and thus 

^(KWk + W u lW\ + W u lWk + W u h\\k) +2rfi||V ft 4||| /! + 2d 2 ||V^||^+2d 3 ||V^I||^ 

< C 2 H\\l h + Ci hl\\% h + C, \\ul\\% h + 2 (77 {ul\ y=e ,ui),ut - ul\ y=e ) % . (50) 

For the last term on the right-hand side of the previous inequality we have 

2{n(ul\ y ^t,u\) % u A h -ul\y =t ) tjfh = 2k(R(ul\ y=( )Q(u' i h ),ul) % -2k(R(ul\y =e )Q(u A h ),ul\y =( )^ 

<0 

< ^{c^eHf^ + C^epy^lf^ + l -\\u 4 h \\l) , 

where we used the assumption (Ai), Young's inequality with e > and the discrete trace inequality (38) 
with the constant C3 > 0. Using the result in (50) we obtain 

l(H\\k + HW* + \\4\\* h + +di\\Vhul\\\ +d % \\V yh ul\\^ h +C4Vy h ulf M , h 

< C 2 \\ul\\\ +C 1 ||« / ;|| 2 ^ +C 5 \\ul\\\ +C 6 \\ui\\\ , (51) 

where C4 := ^3 — kc^C^E can be made positive for e sufficiently small, C5 := C\ + kc^C^e, and Q := kc q j. 
Discarding the terms with discrete gradient, we get 

I (kill* + \\4\\% h + \\4\\% h + \\4\\h h ) < ci(\\ul\\ 2 « k + \\u 2 h \\% h + \\ul\\\ + \\uX h ), 

where C7 :— max{Ci,C2,Cs,C(,}. Applying the Gronwall's lemma to the previous inequality we obtain 

max (|| U j[(0|||, + ||4(0II^ + \\4(f)\\% + \\4(t)\\\) 

< (ll4(0)||| fi + \\4(0)\\% h + \Wl(0)\\% h + \\4{Q)\\l^ T . (52) 
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Finally, integrating (51) over [0,7'] and using (52) gives 

,1||2 , HV7 2-112 , 11*7 ...3||2 



l|V A "il|^ + ||V,^||^ + ||V^||^jd/ 
1 +C 1 Te c 7 T 



< 



(114(0)111, + ll«?(0)||^ + ||^(0)||^ + ||4(0)||| ft ), (53) 



c 8 

where Cg := min{^i , cfe, Q}. The claim of the lemma directly follows. □ 
Lemma 8. Lef {m^,m^,m?,u^} £>e a semi-discrete solution of (2) for some T > 0. 77ien if /ioWs f/ta? 

(ll4(0lll A + ll«£(0lfo + Il4(0ll^) < c, (54) 







max | 

res 

' V ft 4||| + IIV^H 2 ^ + HV^II^ )d/ < C, (55) 



where C is a positive constant independent ofh x , h y . 

Proof. We follow the steps of [19, Theorem 4]. Differentiate (44)-(46) with respect to time, take <p' h 
u' h ,i= 1 , . . . , 3, discard the negative terms on the right-hand side and sum the inequalities to obtain 



1 d 
2dt 



\4\& h + 11411k) +^11^4111+^21^^11^+^1^^411^ 

< Bi M {\ + H) (4,4l>=oL - Bi M \\iil\ y=0 \\\ + (a +« (4,4) 



, 11 — /i i_v — vuwf, 1 v™ 1 r 1 \~ni~~nj jr. 

3 I „4\ .3 1 1 3 „C,.3| ,.4\ .4 .3 1 



As in the proof of Lemma 7, for the first two terms on the right-hand side we have that 

B/ M (i+//)(4,4l v= o) f#;i -B/ M |i4l,=oll| ; ,<c 1 |i4ii| ;i , 

and for the third term 

(a + P)(4,4j, h <C 2 (\\ul\\% h + \\ul\\%J. 

Using the Lipschitz property of 77, together with Schwarz's and Young's inequalities, and assuming the 
structural restriction d, 77 > 0, we obtain for the last term on the right-hand side that 

- (d r T]ul\y = t + d s r]ut,ul\ y= e)y h = -(d r riul\y = i,ul\ y=e ) % - (d s ^utul\ y=t ) % 

< -(d r riul\ y=e ,ul\ y=e ) % +C (^|KII^ + |ll4Ml^) • 
<o 

Choosing e sufficiently small, we get that 

-(d r 7] (4 \y=i , 4)4 \y=t + dsT] (4 !>•=*> 4)4 > 4 \y=i) % < c 3 1 14 1 1 ^ ■ 
Putting the obtained results together we finally obtain that 

^(ll4lll„ + 11411^ + +« , il|V/,4iil ; , +d2\\v yh *l\\h, +*||v >ft 4n^ 

< ^11411^+^(11411^ + n4ii 2 ^) +c 3 ||4lll fc - (56) 
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Gronwall's inequality gives that 

max (|K|||, + \\ul\\% h + \\ul\\% h ) < C 4 (||4(0)||| A + \\u 2 h (0)\\% h + ||^(0)||^). (57) 

In order to estimate the right-hand side in the previous inequality, we evaluate (40)-(42) at t = and test 
with («i(0),M2(0),«3(0))toget 

114(0)111, + Il4(0)||^ + ||4(0)||^ = d!(A ft 4(0),4(0))^ + d 2 (A yh u 2 (0),u 2 (0))# h 
+ d 3 (A yh ul(0),ul(0)).? h — Bi M (Hu\ (0) - ^(0)| y= o,4(0))^ 

+ (au 2 h (0) - J3^(0) , u\ (0) - 4(0))^ . 

Schwarz's inequality and Young's inequality (with e > chosen sufficiently small) together with the regu- 
larity of the initial data yield the estimate 

114(0)11!, + 114(0)11^ + ||4(o)ll^ <c, 

where C does not depend on the spatial step sizes. Returning back to (56), integrating it with respect to t and 
using (57) gives the claim of the lemma. □ 

In the following lemma we derive additional a priori estimates that will finally allow us to pass in the 
limit in the non-linear terms. In order to avoid introducing new grids, grid functions and associated scalar 
products for finite differences in x variable, we will resort to sum notation in this proof. To this end, for 
«/, G J?),, let 8~u.ij, 8yUij, 8~Uij denote the forward and backward difference quotients at xij in x- 

and y-direction, i.e., 

( S x u h)ij 

{5y u h)ij 

Lemma 9 (Improved a priori estimates). Let {uj^u^u^u^} be a semi-discrete solution of (2) for some 
T > 0. Then it holds that 

max (h x h y N f} £ (5+4)2 + hxhy N ^ £ (8+ul) 2 ) < C, (58) 

J£i V i=0 7=0 i=0 ;=0 

T N x -\N y -l j N X -\N V -1 

/ M, £ £ (8 x + 8+u 2 j) 2 dt+ / M, £ E (S/S + 4) 2 df<C, (59) 
70 ;=o i=0 J0 i'=0 7=0 

where C is a positive constant independent ofh x , h y . 

Proof. Following the steps of [19, Theorem 5], introduce a function # G Cq(Q.) such that < # < 1 and let 

#/7 := #b,, G %. Test (17b) with -8~ {& 2 8+ul)ij, (17c) with -S~(t9-?S+4);./, and sum over co h to form 



h x 

u i,j+l — Uij 

hy 



($x u h)ij ■ 

{8yU h )jj : 



Ujj Uj—ij 

h~ x 

Ujj UjJ—1 

h y 
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relations analogous to (45), (46). We get 



i=iy=o i=i 7=o 

= -Bi m £ i (Huj ul )8- (tffi+u&o + <*hy E t i Yji4j8 x -(#i8M)ij 

i=\ i=\ j=0 

•=lj=0 



h I E ^ ij^-^uDij - d 3 Ay E e' rh/}8;4s; (57 Wfe + "*))y 

'=1 7=0 i=l 7=0 

= E )^4v M *) 5 ^ 5 > 3 ^ - aft, EE iyj^M^i^hia 

i=\ i=l 7=0 

+^^EE^^4^(^ 5 , + ^)-/- 

i=l/=0 

Summing the previous two equalities and using the discrete Green's theorem analogous to (34), Schwarz's 
inequality and Young's inequality we obtain 

1 d ' Nl ~ l ^ ^ S N x -lNy-l 

2dt 



^TA h y I El*H + 4l 2 +^ E EiM + 4i 2 )+^.v E E l*W4l 2 

i=0 7=0 i=0 7=0 i'=0 7=0 

N x -lNy-l M x -1 N x -i N y 

+d 3 h y E E {tiS+sf^KBrHd E iM+u/p+cyi, E E(M + 4) 2 

!=0 7=0 i=0 (=0 7=0 

i=0 7=0 i=0 

We rewrite the last term on the right-hand side as 

-k N £\8+(R(ul^M4)M$rs>i.N y ) 

i=0 



* (fiC^W-U) +*(«*+!,* y )5+Q(«f)) (dffi+i^) 

!=0 



i=0 (=0 



< 
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where we used the monotonicity of R and boundedness of Q. To estimate the last term we exploit the 
Lipschitz continuity and boundedness of Q and use the discrete trace theorem so that 

-* £^(^1 «)5+<2(^)(^5+«^) < C 4 \ #?|S+4v/+«?| 

r „iv v -i r. N *-^ N x -lNy 

< I (^5/^,) 2 +g E (d^V) 2 <c 5£ /v E E(M + 4) 2 

Z 1=0 Zfc (=0 i=0 j=0 

i=0 ;=0 /fc i'=0 

Using the latter result in (60), we arrive at 

^fc E E(^4) 2 +^ E L(^4) 2 )+d 2 h y e E (^W« 2 ) 2 

Z Ur 1=0 ;=0 i=0 ;=0 i'=0 ;=0 

+ (rf3-c 5 e)Av e' E ft#8+4) 2 ^ Bi "' HC ' "l! l* 8 ?*} I 2 +c 2 /% e' E (<H + 4) 2 

i=0 y'=0 i=0 i=0 ;=0 

+ (C 3 + qjeJAy e' E (M + »y) 2 + S I W"*) 2 - (61) 

1=0 ;=0 Zfc 1=0 

Applying Gronwall's inequality and integrating with respect to time we obtain the claim of the lemma. □ 

5 Interpolation and compactness 

In this section, we derive sufficient results that enable us to show the convergence of semi -discrete solutions 
of (2). To this end, we firstly introduce extensions of grid functions so that they are defined almost every- 
where in £2 and co and can be studied by the usual techniques of Lebesgue/Sobolev/Bochner spaces. Finally, 
we use the a priori estimates proved in section 4 to show the necessary compactness for the sequences of 
extended grid functions. 

5.1 Interpolation 

In this subsection we introduce extensions of grid functions so that they are defined almost everywhere in £2 
and co. 

Definition 10 (Dual and simplicial grids on £2). Let £2/, be a grid on £2 as defined in Section 3.1. Define the 
dual grid £2^ as 

£2° := {Jf® c £2 | := [*,- h x /2, Xi + h x /2] n£2, x, e £2/,}, 
and the simplicial grid £2^ as 

Slf := {J^ H C £2 | := [ Xi ,x i+1 ]nd, x, e £2,,}. 
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Definition 11 (Dual and simplicial grids on i2 x Y). Let ©/, be a grid on x Y as defined in Section 3.1. 
Define the dual grid (of as 

< := {JSJJJ C £2 X f I ^ := fx, - h x /2, Xi + h x /2] 

x [y ; - - h y /2 , ^ + Ay/2] n A x ? , x, e n A , yj e y A } , 

and the simplicial grid 0)^ as (0® := U ©£, where 

fitf : [(x,-,y ; ), (x ;+1 ,)>;), (**,yy + i)] x HQ x f , i = 0, . . . ,N X - 1 , J = 0, . . . ,N y - 1} , 

< := {3??j | J2£ := [(x i+u y j+1 ), (x i+1 , yj ), (^^i^nflx ?,/ = 0, . . . - 1,; = 0, . . . ,iVy - 1}, 



where [x,y,z] x denotes convex hull of points x,y,z e 



1.2 



Definition 12 (Piecewise constant extension). For a grid function m/, we define its piecewise constant exten- 
sion Uf, as 

■* w= w x€^,^ e ^. (62) 

Definition 13 (Piecewise linear extension). For a grid function Uf, € Sf/, we define its piecewise linear exten- 
sion t?/, as 

u h (x) = Uj + (V h u h ) i+l/2 (x-Xj), x g J^ s , u h e (63) 
while for m/, S i^/, we define it as 

+ Sfuij+ifc+i -x) + (Vyhu h )i + i j+i/ityj -y), x e jsf£. 

The following lemma shows the relation between discrete scalar products of grid functions and scalar 
products of interpolated grid functions in L 2 (Q.) and L 2 (Q. x Y) and follows by a direct calculation. 

Lemma 14. /r /jo/cfo f/iaf 

("/nV/Oz^nxy) = ( u h,Vh),^ h , u h ,v h e ^},, 

5.2 Compactness 

In this subsection we prove our main result. To do this we essentially use the preliminary results shown in the 
previous paragraphs and the results of [13]. Basically, we show the convergence of semi-discrete solutions 
to a weak solution of problem (P). This result is stated in the following theorem. 

Theorem 15. Assume (A1)—(A4) to be fulfilled. Then the semi-discrete solution {uj t , u\, u\, ufy of (2) exists 
on [0, T]for any T > and its interpolate \u\, uj v u\, u\} converge in L 2 (i2), L 2 (£2 x Y), L?(£l x S), L 2 (£2), 
respectively, as \h\ — >• to a weak solution (u\, u 2 , W3, U4) to problem (P) in the sense of Definition 1. 
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Proof. We start off with recovering the initial data. The definition of interpolation of grid functions leads, 

as \h\ ->■ 0, to 

fij(0) — » m? weakly in// 1 ^), 
4(0) -> m° weakly in L 2 (£l;H l (7)), 
4(0) ->■ M3 weakly in L 2 (D.;H l (F)), 
4 (0) M4 weakly in L 2 (£2) . 

Let /i„ be a sequence of spatial space sizes such that \h\ — > as « — > 00, Consequently, we obtain a 
sequence of solutions {kJ , , 4 ; M /!„ } of ( 17) defined on the whole time interval S. 

Let us pass to the limit \h\ — > in the ODE. Note that f]{u\ \y=e,ul ) — q weakly in L 2 (S;L 2 (£l)), and q 
still needs to be identified. The way we pass to the limit in the ODE is based on the following monotonicity- 
type argument (see [21]): using the monotonicity of 77 w.r.t. both variables, we can show that u\ is a Cauchy 
sequence, and therefore, it is strongly convergent to u 4 . 

Now, it only remains to pass to the limit in the PDEs. Note that the weak formulation contains a nonlinear 
boundary term involving rj(-, •). Exploiting the properties of the interpolations of grid functions we deduce 
that the same a priori estimates hold also for the interpolated solution (see also [13]). On this way, we obtain 

{u\ n } is bounded in L°°(0, T;L 2 (Q.)), 
{u\ n } is bounded in L 2 (0, T;H X (£2)), 
{alj is bounded in L°°(0, T;L 2 (Q.)), 
{a 3 hn } is bounded in L°°(0, T;L 2 (Q.)), 
{u 4 h J is bounded in L°°(0, T;L 2 (Q.)). 

Hence, there exists a subsequence of h„ (denoted again by /?„), such that 

a l hn u l weakly in L 2 (S;H l (£2)), 
u\ n u 2 weakly in L 2 (S;L 2 (£2)), 
u\ n m 3 weakly in L 2 (S;L 2 (£2)), 
u\ n -± u A weakly in L 2 (S;L 2 (£2)). 

Since 

ll"ft„llL2(s,Hi(a)) + Il^"/,jz2(s,z.2(n)) < c > 
Lions-Aubin's compactness theorem, see [14, Theorem 1], implies that there exists a subset (again denoted 
by u\ ) such that 

u\ — > it 1 strongly in L 2 (5x£2). 

To get the desired strong convergence for the cell solutions k 2 , , v? h , we need the higher regularity with 
respect to the variable x, proved in Lemma 9. We remark that the two-scale regularity estimates imply that 

W^ljl^iS-H 1 {n.fli (Y))) + \\ul n \\L 2 {S;H l (n.H i {Y))) - C - 

Moreover, from Lemma 8, we have that 

\\dtU 2 hn \\ L 2 {Sx£lxY) + \\dtul n \\ L 2(sxnxY) ^ c - 
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Since the embedding 

H l (Q.,H l (Y))^L 2 {Q.,HP{Y)) 

is compact for all A < j3 < 1, it follows again from Lions- Aubin's compactness theorem that there exist 
subsequences (again denoted uj^ , u\ ), such that 

(fij^)— >(hV) strongly in L 2 (5xL 2 (fi,#(y)), (65) 
for all ^ < j3 < 1 . Now, (65) together with the continuity of the trace operator 

H p {Y)^L 2 (dY), for I<j3<l, 
yield the strong convergence of , w| until the boundary y = 0. □ 

6 Numerical illustration of the two-scale FD scheme 

We close the paper with illustrating the behavior of the main chemical species driving the whole corrosion 
process, namely of H2S(g), and also the one of the corrosion product - the gypsum. To do these computa- 
tions we use the reference parameters reported in [5]. 




2 4 6 8 10 2 4 ,, 6 8 10 2 4 . 6 8 10 2 4 ,, 6 8 10 



Figure 1 : Illustration of concentration profiles for the macroscopic concentration of gaseous H2S (left) and 
of gypsum (right). Graphs plotted at times t G {0,80, 160,240,320,400} in a left-to-right and top-to-bottom 
order. 

Figure 1 shows the evolution of u\ (x,t ) and u^x^t) as time elapses. Interestingly, although the behavior 
of «i is as expected (i.e., purely diffusive), we notice that a macroscopic gypsum layer (region where M4 is 
produced) is formed (after a transient time t* > 80) and grows in time. The figure clearly indicates that there 
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are two distinct regions separated by a slowly moving intermediate layer: the left region - the place where 
the gypsum production reached saturation (a threshold), and the right region - the place of the ongoing 
sulfatation reaction (Id) (the gypsum production has not yet reached here the natural threshold). The precise 
position of the separating layer is a priori unknown. To capture it simultaneously with the computation of 
the concentration profile would require a moving -boundary formulation similar to the one reported in [4]. 
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